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ABSTRACT 



A discrete potential element approach to subsonic numer- 
ical lifting surface theory has been developed and shown to 
be practical in predicting the nonsteady loading on 
harmonically oscillating, medium to low aspect ratio wings. 

A unique method of including the wake effect in the wing 
kernel function matrix prior to solution of the singular 
integral downwash equation was devised, thus greatly 
simplifying the velocity potential formulation. In addi- 
tion, termination of the effective wake a finite distance 
downstream of the wing was investigated, with wing loading 
found to converge to within one per cent in an effective 

4 

wake length of four root chords. 

This discrete element method has also been extended to 
the case of an oscillating wing, cantilevered from a cylin- 
drical fuselage, to investigate nonplanar interference 
effects. This interference in wing loading, while of rela- 
tively small magnitude, does exist in both pressure amplitude 
and phase angle distributions, and is, therefore, of impor- 
tance in three-dimensional stability analyses of wing/body 
configurations . 
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I. INTRODUCTION 



The aerodynamic analysis of the forces over an oscillat- 
ing wing, situated in a steady flow, is the basic problem 
considered by nonsteady lifting surface theory. The deter- 
mination of these nonsteady pressure distributions is of 
prime importance in the consideration of stability. A basic 
objective of flutter analysis, for example, is to determine 
the aerodynamic loading on a wing oscillating with small 
amplitudes in a definite mode and with a given frequency. 

Thus the problem is formulated on the basis that the motion 
and deformation of the lifting surface are either known or 
composed of known elemental modes, and is amenable to har- 

4 

monic modal analysis. A downwash integral equation, based 
on the transformed wave equation for linearized potential 
flow, relates loading (either in terms of velocity potential 
or pressure) over the oscillating wing to the normal velocity 
induced by the wing motion, and must be solved numerically 
to obtain this wing loading. 

The analysis of this paper will deal primarily with the 
subsonic case, this being the one in which the largest 
number of applications of lifting surface theory to engineer- 
ing problems are to be found [1] . Most of the work in sub- 
sonic nonsteady lifting surface theory is based on the 
development of the general kernel function approach by 
Kussner in 1940 [2], which was first formulated successfully 



for the computer by Watkins, Runyan, and Woolston in 1955 
[ 31 . More recent techniques have been developed by Stark in 
1964 [4] and Laschka in 1963 [5] . Solution of the singular 
downwash integral equation by the usual lifting surface 
methods is accomplished through assuming the loading to be 
a series of preselected functions with unknown coefficients, 
and then determining these by satisfying the normal velocity 
condition with some sort of collocation technique. The 
downwash and the loading are related by kernel functions, 
which are themselves singular, and which must be numerically 
evaluated. 

In spite of the fact that lifting surface theory has 
enjoyed a great deal of success, certain shortcomings, 
discussed in Refs. 1, 6, and 7 , do exist: 

(i) Lifting surface theory is not as computationally 
economic as steady or two-dimensional methods. 

(ii) Wing loading functions, which incorporate the 
proper singular behavior, have to be assumed a priori. 

(iii) Uncertainty exists as to the method of locating 
collocation points. 

(iv) Techniques for analyzing control surface effects 
have not been adequately developed, partly because of (ii). 

(v) Lifting surface theory has been primarily restricted 
to planar surfaces. 

In order to attack these problems, Houbolt [6] proposed that 
the wing loading be represented by concentrated pressure 
loads in the form of Dirac delta functions. Using this 
loading representation, the kernel function and downwash 
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equation were formulated and analyzed for certain restricted 
cases. However, no general numerical solution of the lifting 
surface problem was developed. 

The present study is concerned with developing such a 
method of solution to the subsonic, nonsteady lifting surface 
problem, in terms of discrete loading elements. However in 
this case wing loading is expressed in terms of velocity 
potential, rather than pressure, since this formulation 
allows a more direct approach to nonplanar configurations. 
Haviland [8] reports using a similar approach to the steady 
planar case with consistent results. 

The aim of this discrete potential element method is to 
provide a new approach to nonsteady, subsonic lifting 
surface theory which eliminates, or minimizes, the problems 
previously indicated. Theoretical development of this 
approach was made from basic principles, and the resulting 
formulation applied in a computer program which analyzes 
harmonically oscillating wings in subsonic flow. A unique 
way of including the wake effect was also developed which 
greatly simplified the velocity potential solution of the 
numerical problem. As an extension of the theory, and an 
indication of the versatility of this discrete element 
approach, a second computer program was developed to analyze 
interference effects on oscillating wings mid-mounted from 
steady bodies. To the author’s knowledge, this type of 
wing/body interference investigation has not previously been 
presented for the nonsteady case. 
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II. GENERAL PROBLEM DEVELOPMENT 



Following the method of development by Garrick [9], the 
governing equations and boundary conditions are linearized 
by application of small perturbation theory. 



A. POTENTIAL FLOW EQUATIONS 

The governing equations for adiabatic, irrotational, 
inviscid flow, as developed in Ref. 10 are: 

Continuity 

^-+V-pti = 0 (2.1) 

3 1 



Momentum 

3U ->■ -* 1 

“ + u • VU = - - VP 

a t p 



( 2 . 2 ) 



Energy 

p 

— = constant (y the specific heat ratio) (2.3) 

P Y 

Since the flow is irrotational, the curl of the velocity 
vector vanishes (Vxi5=0) and a velocity potential function 
can be defined, such that 

tl = V$ (2.4) 



If further, a perfect gas is assumed, then from the energy 
equation the well known relationship 

3P/3 p = a 2 (2:5) 

is obtained for constant entropy flows. 
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The continuity and momentum equations can now be rewritten as 



i £!£. + 

p Dt 



ii + HI 
9t 2 



= 0 



Vp 



where ^ is the substantial derivative defined by 



9t 



( 2 . 6 ) 

(2.7) 

+ U • V| 



Combining equations (2.6) and (2.7) , a single equation 
in $ is obtained. 



V 2 $ - -4 



at 



+ U • V 



= 0 



or 



i B?* . o 

a Dt 



( 2 . 8 ) 



This is the governing partial differential -equation for 
general nonsteady, nonviscous, potential flow of a perfect 
gas. Equation (2.8) is not limited to small disturbances, 
and has the form of the classic wave equation in terms of 
the substantial, or convective derivative. Thus, the flow 
disturbance represented by the velocity potential is con- 
vected by the local fluid, or stream velocity, and propogated 
as a wave which spreads at a rate equal to the local speed 
of sound. 

In this general form, equation (2.8) is highly nonlinear 
due both to the quadratic terms in the convective operator 
TJ*V, and to the interdependence of the velocity potential 
and the local speed of sound. Boundary conditions, in 
general, will depend on the location and form of the moving 



18 



body, on discontinuities along streamlines, on shock waves, 
and on flow conditions at infinity. There exist no general 
methods for obtaining solutions to equation (2.8), and it 
is necessary to consider small perturbations, to linearize, 
to reduce the number of equations, or to try schemes of 
successive approximations to attack the general problem. 

B. SMALL PERTURBATION THEORY 

Most of the theoretical developments in aerodynamics, 
including nonsteady aerodynamics, are based on the concept 
of small perturbations. There are two aspects to the prob- 
lem of flow past a body creating small disturbances in an 
undisturbed mainstream, linearization of the governing 
differential equations and further linearization of the 
boundary conditions. 4 

Linearization of the governing equations requires the 
problem to be posed with a compressible fluid in a uniform 
stream flowing with velocity in the positive x direction. 
Small disturbances in the flow are defined by 

5 = (U^+u) i + v j + wk = U^i + u 

P=P 0O +P $ = ^ $ 

P = p^ + p with u = V<5> 

The perturbation velocities are considered small with respect 
to U , a , and U -a , thus ruling out transonic flow. The 
perturbation density and pressure are also restricted, so 
that 
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P/Poo and P/Poo <<] - 



Considering <j> as the perturbation velocity potential, 
the governing equation (2.8) becomes 

,2 



u 2 I 1 
V <J> - — 

a„ 



A + u — 
at °° ax. 



$ = o 



(2.9) 



where higher order terms in the U'V operator have been 
neglected. This linearized governing equation may now be 
applied to flow problems within the initial small perturba- 
tion assumptions. Although equation (2.9) was developed for 
a uniform flow in the positive x direction relative to a 
space-fixed coordinate system, it also applies to the case 
of a wing moving with uniform velocity in the negative x 
direction with the coordinate system fixed to the wing. 

4 

Within the concept of small perturbations, the momentum 
equation becomes 



5u 

3t 



+ 



u |5 

00 3x 




( 2 . 10 ) 



Substituting the perturbation velocity potential (u=V<j>) 
into this equation, a linear relationship between potential 
and pressure is obtained. 



P = -P, 



— + u — 

at ~ 9x 



9 



(2.11) 



Therefore p satisfies the same differential equation (2.9) 
as does <j>. 

If harmonic motion were considered, the perturbation 
variables would have the form 
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t i / \ loot — / \ lwt 

<f> = <j>(x,y,z)e , p = p(x,y,z)e 

u = u(x,y , z)e^ wt , etc. 

and the governing equation would become 

V 2 4> _ A ( Uoo + ±oj) 2 4> = 0 (2.12) 

a oo 

with the relationship between pressure and velocity potential 



p = -p (U 

r ^oo < 




ico ) <p 



( 2 . 13 ) 



In the above equations the time dependence e ia)t has been 
factored out. This relationship between p and <{> can be 
integrated to give 



4>(x,y >z) 



-iwx/U 

CO 

e 



P U 

~ CO CO 



x iwC/U^ 

/ p(€>y>z)e °° d£ 



( 2 . 1 - 4 ) 



where far ahead of the disturbance, the perturbation poten- 
tial is assumed to be zero. 



C. LINEARIZED BOUNDARY CONDITIONS 

The small disturbance assumption, resulting in the 
linearized governing differential equation (2.9), implies 
small deviations from the uniform flow. Therefore it is 
also necessary to linearize the boundary conditions for the 
physical problems considered, such as properly oriented thin 
wings and slender bodies. These boundary conditions consist 
of the mathematical formulation of several physical and phe- 
nomenological statements which are normally grouped as 
follows: 
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(i) Surface boundary conditions. No flow can pass 
through the wing or body surface, that is to say the total 
flow is tangential to the body surface. 

(ii) Edge conditions. Sufficient viscosity is present 
in the "nonviscous" flow to determine the flow pattern near 
sharp edges. Thus in subsonic flow the well known Kutta 
Condition can be applied to the wing trailing edge, requiring 
that the surface pressure difference remains continuous 
near, and vanishes at, the trailing edge. A similar condi- 
tion should be satisfied at the side edges. Through the 
method of matched asymptotic expansions, Landahl [11] veri- 
fied that at trailing and side edges the pressure difference 
will approach zero with infinite slope due to weak singu- 
larities in the first derivative. Conditions at the leading 

* 

edge have not been fully explored, but for small disturbances 
to rounded leading edges it is sufficient to require that 
the total integrated force be finite and the singularities 
in the pressure distribution be of the proper order [91 • 

This latter requirement is fulfilled with the pressure 
varying as the inverse square root of the distance downstream 
of the leading edge [12] . 

(iii) Wake conditions. The free vorticity shed from 
the trailing edge in the wake is such that its circulation 
together with the bound circulation vanishes in accordance 
with the Helmholt z-Kelvin Theorem. It is assumed that the 
shed wake remains where it is formed, floats without mutual 
interference along streamlines, and forms a continuous sheet 
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of discontinuity coplanar with the wing in the direction of 
flight. Edge effects and roll-up of the sheet at infinity 
are ignored. 

(iv) Conditions at infinity. The flow is considered 
uniform at infinity, and perturbation waves are required to 
propagate away from the sources of disturbances and to 
behave properly at infinity. 

In addition to the above, for supersonic flow account 
must be taken of zones of influence and dependence, as in 
the method of characteristics. In supersonic flow there are 
also other problems which do not arise in the subsonic 
analysis, such as the difference between subsonic and super- 
sonic edges. 

To formulate the main boundary conditions, a thin wing 
is considered lying in the x-y plane, creating small dis- 
turbances in a uniform stream flowing in the positive x 
direction. The linearized tangential flow condition on the 
wing surface z g takes the form 



The normal fluid velocity w can be expanded in a power series 
about z = 0, so that 



w(x,y,z g ,t) = + U 



00 3x 



(2.15) 



w(x,y,z g ,t) =w(x,y,0,t) + 




z=0 



z 



s 




3 2 w| -2 . 

— z + 



2 s 

3z /z=0 S 



( 2 . 16 ) 
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In addition, problems of camber and thickness can be separ- 
ated from wing motion and angle of attack in the linearized 
formulation. Therefore, neglecting higher order terms in 
equation (2.16), the tangential flow. condition becomes 

3z 3z 

w(x,y,0,t) = ^ + U ro (2.17) 

where z (x,y,0,t) represents the mean camber surface posi- 
tion of the wing. ' 

Since w is an even function of z (its value is unchanged 
for the top or the bottom side of the mean surface), the 
velocity potential <j) is an odd function of z, as is the 
perturbation pressure p through the linear relationship of 
equation (2.11). Therefore, in the plane of the wing 6 
equals zero outside of the wing and wake. 4 > does not equal 
zero in the wake because of the discontinuity in the u 
component of the perturbation velocity across the wake. The 
perturbation pressure is also zero in the x-y plane, except 
at the lifting surface where the pressure difference Ap 
between the top and bottom surfaces is given by 



Ap(x,y,t) = p(x,y,0-,t) 
= 2p(x,y ,0,t) 



p(x,y,0+,t) 



“2p„ 



_ 3 _ 

3t 



+ U 



_ 3 _ 

3x 



4 > 



( 2 . 18 ) 



Figure 1 summarizes the formulation of these boundary condi- 
tions for a lifting surface in the x-y plane. 

These linearized boundary conditions along with the 
linearized governing differential equation (2.9) can be used 
to attack a wide range of perturbation problems in subsonic 



2b 
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and supersonic flow. Miles [131 has shown that a sufficient 
condition for linearization of the equations for the pertur- 
bation potential is any one or combination of the factors 



|m 2 - 1 | » 0 (6 2/3 ) 

k » 0 (6 2/3 ) 



^»0 ( 6 1/3 ) 



where 6, M 6, k6 , and kM 6 must be sufficiently small. 

J CO J J CO v 

The perturbation equation is nonlinear only when the 
following conditions are jointly satisfied 



M 2 - l| = 0 (6 2/3 ) 
k = 0 (6 2/3 ) 



1 _ n ( 

AR - 0 (6 } 



6 is here intended to be a thickness, angle of attack, camber, 
or motion parameter serving as a measure of the flow 
disturb ance . 
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III. LIFTING SURFACE THEORY 



A. BASIC APPROACH 

The basic nonsteady problem considered in lifting surface 
theory is that of a linearized thin wing harmonically 
oscillating in the normal direction in some prescribed 
deflection mode, creating small disturbances in the uniform 
mainstream. Particular solutions to the linearized governing 
differential equation (2.9) are required which give pertur- 
bation velocity potential or perturbation pressure distri- 
butions on the wing, and which satisfy the essential boundary 
conditions discussed in Section III-C. Such a direct 
solution has been possible for only a few special cases, 
such as two-dimensional incompressible flow as in Refs. 14 
and 15. For the more general three-dimensional case, an 
integral equation formulation must be employed which makes 
use of a downwash kernel function, analogous to an influence 
coefficient development. 

Considering the time dependence e la> in the harmonic 
analysis as being factored out, the integral equation relating 
the downwash amplitude to the wing loading has the form 

w(x,y,0) = // L(5,n,0)K(x-£, y-n,0)d£dn (3-D 

S 

The loading L can either be the velocity potential or the 
pressure distribution over the surface, and K is the kernel 
function which denotes the normal velocity (downwash) on the 
wing at (x,y,0) due to a unit acoustic singularity located 
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at (£,n,0). In other words, the surface is modeled by a 
distribution of radiators, formulated in terms of either 
potential or pressure which satisfy the linearized governing 
equation 



V 2 4> “ (U*, 37 + lw ) 2 * = 0 (2.12) 

a <» 

w is the downwash at point (x,y,0) on the wing, determined 
by the harmonic motion of the surface in fulfilling the 
tangential flow requirement. All of these quantities are 
complex in the spatial domain due to the nonsteady nature of 
the problem. 

Employing the velocity potential approach, the kernel 
function is defined by 

K = w(x ,y , 0 ) = d<p/dz 4 (3.2) 



where 0 is the solution to equation (2.12) for an acoustic 
radiator of unit strength. In this formulation, the inte- 
gration (3-1) must be carried out over the surface of both 
wing and wake, since the potential does not go to zero in 
the wake region. This complicates the integration process 
and has proven a major disadvantage of the potential approach. 

If the wing is modeled by a pressure distribution, the 
kernel function is denoted by 



K = w(x,y ,0) 
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from equation (2.1*0 where p is now the solution to equation 
(2.12) for an acoustic radiator of unit strength. The 
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integration here need only be carried out over the wing 
surface, since the pressure goes to zero in the wake. 

Because of this, methods based on the pressure formulation 
have been more widely employed. This kernel function (3-3) 
is more complicated than that developed from equation (3.2) 
and has not been integrated in closed form due to singulari- 
ties at (y-n) = 0 and (x-£)>. 0 . However these singularities 
have been isolated and expressed in forms which can be 
handled by numerical procedures [3] . 

The acoustic singularities used to model the wing surface 
in subsonic flow are harmonically pulsating doublets (dipoles) 
oriented in the z, or wing-normal, direction. Doublets are 
employed in order to represent the pressure difference 
between the two surfaces of the wing, while in supersonic 

4 

flow acoustic sources, or monopoles, are used due to the 
independence of the two wing surfaces in that flow regime. 
Kiissner [2] developed the general formulation of the pressure 
doublet kernel function, through use of the Lorentz trans- 
formation, in order to apply solutions to the classical wave 
equation of physics 



a 3t 

to the governing differential equation of lifting surface 
theory (2.9). This method of attack is used in Appendix A 
to develop the kernel functions used in this report. 
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B. NUMERICAL SOLUTIONS TO LIFTING SURFACE THEORY 



To find the potential or pressure loading on the lifting 
surface actually requires the inverse solution of equation 
(3-D) since the downwash velocity distribution is deter- 
mined by the type of surface motion prescribed through 
equation (2.15). The numerical methods that have been 
developed in subsonic lifting surface theory rest on this 
inverse solution to the singular integral equation. The 
usual procedure is to assume the loading to be a series of 
preselected functions with unknown coefficients and then 
to determine these by satisfying the downwash velocity 
distribution exactly in a set of collocation points [16, 17], 
or approximately in the least squares sense in a larger set 
of control points [4, 18] , or by satisfying certain integral 
relations derived from variational procedures [19, 20]. 

Ashley and Landahl [21], and Landahl and Stark [1] have 
presented survey papers summarizing the status of numerical 
lifting surface theory, and presenting formulations of the 
integral downwash equation and the kernel function in the 
various methods of attack, such as: 

(i) Velocity potential. As previously discussed, the 
kernel function used in this formulation, while singular, is 
much simpler than that employed with the pressure loading 
approach. However, integration must be carried out over 
both wing and wake. The velocity potential formulation has 
been employed by Jones [22] for = 0 in 1952. 

(ii) Acceleration potential. Formulating this approach 
in terms of an acceleration potential defined by c = -2 ip [231, 

ir 
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allows the use of the same kernel function as in the velocity 
potential formulation, but removes the necessity of integrat- 
ing over the wake, since c^ = 0 there. However, the solution 
is not unique, since multiples of eigen solutions with 
1^- = 0 on the wing may be added, and integration must be 

o Z 

extended into the region ahead of the leading edge to achieve 
uniqueness. This formulation has been used for two- 
dimensional cases and also for certain three-dimensional 
cases [24] . 

(iii) Pressure formulation. The original development 
by Kiissner [2] has been formulated for the computer by 
Watkins, et al [3] and by Richardson [ 1 6 ] . This procedure 
provides direct determination of the pressure with integra- 
tion required only over the wing surface, but the kernel 
function is highly singular and needs to be evaluated through 
numerical quadrature, thereby increasing computer time 
considerably . 

Various refinements of the basic approaches discussed 
above have been developed, such as the integrated acceler- 
ation potential by Stark [4] and the advanced velocity 
potential [1] , but these do not change the basic problem 
formulation, nor the basic advantages and disadvantages of 
the continuous loading function methods. 

In lifting surface theory the loading is expressed as 
a linear combination of preselected functions, with coeffi- 
cients to be determined by satisfying the tangential flow 
condition through the integral equation (4.1). In selecting 
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the loading functions, known results from two-dimensional 
incompressible flow theory have been used. The loading is 
separated into chordwise and spanwise functions, so that a 
pressure function element would be represented as 

Pmn = f m < 5 ) e n (n) 

These loading functions have been represented by' Fourier 
expansions [17] , power series [25> 26] , Tschebycheff 
polynomials [27, 28], Legendre polynomials [^], and so on. 

The basic requirement is that these functions themselves 
satisfy the applicable boundary conditions from Section II, 
such as leading, side, and trailing edge behavior. 

C. SOLUTIONS BASED ON DISCRETE LOADING LINES 

In steady flow, the lifting line theory* of Prandtl 
employed a discrete vortex lifting line, with variable 
strength, placed at the quarter chord line. Wieghardt [29] 
extended this method for rectangular wings by using several 
lifting lines, while Rubbert [30], Dulmovits [31] 5 and 
Hedman [32] assumed constant strengths of the lifting lines 
in subintervals in a vortex lattice approach. 

The vortex lifting line in steady flow corresponds in 
the nonsteady case to a potential doublet strip with strength 
varying harmonically in the streamwise direction. The use 
of a discrete lifting line approach to the nonsteady case 
was first proposed by Jones [22] for the incompressible 
case, and by Runyan and Woolston [25] for general subsonic 
flow. 
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As in the steady case, this lifting line approach for 
nonsteady motion has been extended to the doublet lattice 
method by Stark [1], and by Albano and Rodden [331. The 
downwash from the doublet lattice strip is closely related 
to the kernel function discussed previously in this section. 
Since the surface loading is replaced by discrete doublet 
strips, it is not necessary to assume the loading functions, 
as required in the continuous function approach to lifting 
surface theory, and the integral downwash equation (3.1) is 
replaced by a matrix equation relating downwash on the wing 
to doublet lattice strength. This approach then removes 
one of the shortcomings to general lifting surface theory, 
but still retains the remaining features. 

D. NONPLANAR CONFIGURATIONS 

Limited work has been done in the area of general non- 
planar configurations, since interference investigations 
have been primarily involved with combinations of lifting 
planar surfaces. The kernel functions for interfering planar 
surfaces, as developed by Davies [3^1 and by Landahl [35] > 
become more complicated but do not contain any new singulari- 
ties. In the selection of suitable loading functions, the 
behavior at the intersection of lifting surfaces must be 
taken into account. The T-tail has been treated by Stark 
[4] and by Davies [35] > while calculations for a delta wing 
with folded tips have been given by Vivian and Andrew [36]. 
Lashka [5] has analyzed the effects of wing tip pylons as 
well as interfering planar surfaces [37] , while Rodden [38] 
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has included wing/body and wing/pylon effects in nonsteady 
wing loading in subsonic flow. 



IV. DISCRETE POTENTIAL ELEMENT DEVELOPMENT 



A. WING ANALYSIS 

1 . Basic Formulation 

The integral equation relating downwash to wing 
loading (either velocity potential or pressure) in lifting 
surface theory, as discussed in Section III, is 

w(x,y,0) = // L(£ ,n,0)K(x-£,y-n,0)d£dn (3-D 

S 

where harmonic motion is assumed so that 

w(x,y,z,t) = w(x,y,z)e iwt 
L ( x ,y , z , t ) = L(x,y ,z)e lwt 

4 

In the subsequent development, the time dependence will be 
considered factored out, so that all symbols will refer to 
complex valued spatial variables. The downwash w is deter- 
mined by the boundary condition of no flow through the wing 
surface as given for general motion by equation (2.7) , and 
here for harmonic motion 

9 z 

w = U — - + iwz (4.1) 

00 dX S 

where z g is the amplitude of the surface motion prescribed 
for the problem. The form of acoustic radiator used to 
model the wing is taken as a potential doublet, with axis in 
the z, or wing normal, direction. This dipole singularity 
is employed because of the requirement to sustain a pressure 



35 



difference across the wing in subsonic flow. K is the kernel 
function relating the downwash at point (x,y,0) on the wing 
due to a unit potential doublet located at (£,ti,0). 

In normal lifting surface theory, the loading L would be 
represented by continuous chordwise and spanwise functions 
which satisfy the surface boundary conditions developed in 
Section II. However, in the discrete potential element 
approach this loading is represented by a network of Dirac 
delta functions formulated in terms of the perturbation 
velocity potential. Thus the velocity potential distribution 
over the wing is replaced by a series of point functions in 
the form of harmonically oscillating potential doublets with 
axes in the z direction. In addition, the wake must also 
contain a network of these doublets since the velocity 

4 

potential is not zero there (Section II) and the integration 
region of equation (3-1) must include both wing and wake. A 
representation of this potential doublet grid is shown in 
Figure 2. 

In the discrete formulation, the integral equation (3*1) 
is replaced by the matrix equation 

{w} = [K] {$} (4.2) 

t h 

where w^ is the downwash at the l point on the wing, and 
4> • is the strength of the j potential doublet on the wing 

J 

and wake in the grid of Figure 2. The location of the down- 
wash control points are coincident with the potential doublet 
positions at the center of each wing control box. The choice 
of this location will be discussed in Section IV. A. 4. In- 
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modeling the wing, the velocity potential and the pressure 
are considered constant over each of the control boxes, with 
approximations to continuous distributions given by these 
values at the control box midpoints. 

The kernel function matrix elements represent the down- 
wash due to potential doublets of unit strength which satisfy 
the governing linearized differential equation 

V 2 <f> " (U M 3^ + iw)2( ^ = 0 (2.12) 

a TO 

As developed in Appendix A, this kernel function has the 
form 

K., = (1+i — ^-R)exp{i — ^ [M w (x -x,)-R]}(4.3) 

13 HttR 3 a B 2 a 1 J 

OO OO 

where 

In this formulation the kernel function is singular only 
when the downwash control point is coincident with the poten- 
tial doublet location (i=j). The handling of this singular- 
ity is vital to this development and is discussed in detail 
in Section IV. A. 4. 

The discrete loading element approach to the nonsteady 
wing problem requires the solution of the system of linear 
complex equations (4.2) to determine the velocity potential 
vector {$). The downwash vector {w} is specified by the 
wing motion through equation (4.1), so that the solution has 
the form 
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(4. 4) 



{*} = [K] _1 {w} 

where { <f> } and {w} cover both wing and wake regions. The 
resultant velocity potential distribution is used to deter- 
mine the pressure loading over the wing through the rela- 
tionship 

Ap = -2p oo (U oo + itoH (4.5) 

taken for the harmonic case from equation (2.18). 

2 . Wake Effect 

One of the historic drawbacks to the potential 
approach to lifting surface theory has been the necessity 
of including the wake region in the solution of the integral 
equation (3.1). This equally complicates the discrete 
element approach since the matrix equat ion «( 4 . 2 ) is required 
to include the effect of the wake potential dipole grid. 

Even though the effective wake is terminated at some repre- 
sentative length behind the wing, as is commonly done in 
steady flow problems, the size of the resulting potential 
strength vector and kernel function matrix would severely 
restrict, or entirely preclude, the use of this method even 
on large modern computers. 

Following a suggestion by Professor R. E. Ball of the 
Naval Postgraduate School, a careful examination of the 
boundary conditions governing harmonic wing motion, as 
depicted in Figure 3, was made. In the wake region the 
pressure is zero, but the velocity potential is not zero 
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because of the discontinuity in the u component of the per- 
turbation velocity. <J> and p are related by 



p = - p « (U ~ dJ + ito) * 



(2.13) 



In the wake, therefore, this relationship becomes 



If + I-O- = 0 



(4.5) 



or 



21 = _ i 



w 



3x 



U 



which can be integrated to give 



d> = <J> e 
Y w Y m 



. to / s 

-irr— (x -X ) 

U w m 

oo 



(4.6) 



<f> is the strength of the wake velocity potential at 
( x w*^i*0)» and 4> m is strength of the velocity potential 

at the boundary between the wing and the wake ( x m >yj_>°) 
along the line of integration y=y^. In the discrete poten- 
tial element formulation, <j> becomes the strength of the 
last wing dipole on the appropriate chord line, with (x -x m ) 
the distance between this dipole and the wake dipole , 
as is shown schematically in Figure 4. 

The wake potential strength distribution is therefore 
expressible as a function of the wing potential distribution 
and can be included in the wing doublet kernel functions. 
This relationship can be represented by 

(<P W ) = [A] {cf>} 
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where {<t»} is the wing velocity potential vector and the 
appropriate elements of [A] express the relationship of 
equation (4.6) with respect to the wake potential vector 

{(b }. The downwash matrix equation (4.2) therefore becomes 

w 

{w} = [K 1 ] (4> > + [K 2 ] 

= [K 1 ] {cf>} + [K 2 ] [A] (<t>> 

= [k 1 +k 2 a] {<t>} 

The size of this matrix equation has thus been limited to 
that necessary to represent only the wing grid, making the 
discrete potential element approach much more tractable on 
the computer. 

3 • Matrix Equation for Symmetric Motion 

In considering a wing undergoing harmonic oscilla- 
tions which are symmetric with respect to the x-z plane, the 
wing/wake planform can be represented as shown in Figure 5* 
The x-z plane, being a plane of symmetry, represents a no- 
flow surface, or reflection plane, for the disturbances 
caused by the symmetric motion of the full-span wing. View- 
ing the physical problem in a different, but equally valid 
light, the half-span wingQ, with its root along the x axis, 
can be considered undergoing oscillations in the presence of 
an infinite wall coincident with the x-z plane. Therefore, 
the motion of the wing need only be prescribed for wing Q, 
while the function of the virtual wing 0 and wake © is to 
establish the no-flow 'condition at the x-z plane. 

The full wing/wake system is represented by the matrix 
equation 
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where subscripts denote the areas of Figure 5 and . is the 
kernel function matrix for downwash at points in area (3) due 
to doublets located in area Q). For symmetric motion 



= and { 4 ^} = { 4> 2 } 



so the required downwash equation is reduced to 



{w-j^} 



[Kn + Ki3 ; 



K 12 +K ld 




\ ) 

4 

In the previous section, the relationship between wing and 
wake potential distributions was developed, and can here be 
expressed as 



{$ 2 } = [A] 



Therefore, the final wing downwash matrix equation has the 
form of equation (4.2) 

{w 1 } = [K] (4.7) 



but the kernel function matrix is here formed from the 
following expression 

[K] = [K n ] + [K 13 ] + [K 12 +K lij ] [A] (4.8) 

Thus, through use of motion symmetry and the boundary condi- 
tion in the wake, the effective integration of equation (4.7) 



45 



need only be carried out over the half-span wing surface Q. 
Acknowledgement of the full physical problem Is obtained 
through formation of the kernel function matrix via equation 
(4.8), which incorporates effects of the potential dipole 
grids in all four areas of the wing and wake depicted in 
Figure 5- 

It should be noted that antisymmetric motion of the wing 
about the x-z plane requires that 

{<£ 3 } = and {<f>^} = -{<f> 2 } 

The final kernel function matrix is therefore formed by 

[K] = [K 1X ] - [K 13 ] + [K 12 -K l2| ] [A] 

Since any general harmonic motion of the wing can be 
expressed as a combination of symmetric and* antisymmetric 
modes, this approach has general application and need not 
be restricted to the symmetric case considered in detail 
here . 

4 . Doublet Singularity 

The location of the collocation points required to 
satisfy the no-flow condition in lifting surface theory is 
historically based on two-dimensional steady flow theory, 
Vortex lattice methods such as Rodden's [ 38 ] arrange the 
vortex strip on the local quarter chord of the control box, 
with the downwash collocation point centered on the local 
three-quarter chord line as in two-dimensional thin airfoil 
theory. Houbolt [6] proposes using this local quarter chord, 
three-quarter chord control box grid in conjunction with 
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concentrated pressure loads. This type of grid network was 
tried by the author in an early form of the wing analysis 
computer program with unsatisfactory results. The approach 
was found not to converge, but to be very sensitive to grid 
size; that is, to the distance between the potential dipole 
and its associated control point. This is much like the 
sensitivity that a continuous loading method experiences 
when a control point is located too close to a wing edge 
111 ] . 

The control or collocation points were subsequently 
placed at the center of the wing control boxes (Figure 2) 
coincident with the potential doublet locations. In this 
way, the above mentioned sensitivity to grid spacing was 
removed, but the value of the upwash of a doublet at its 

4 

own control point had to be determined. As can be seen 

from equation (4.3), the kernel function (K. .) is singular 

J 

at the doublet location (i=j). This singularity in the 
upwash from a doublet is pictured in Figure 6a, where the 
doublet produces an infinite upwash at its location, but 
finite and decreasing downwash in the plane of the wing as 
the distance from the doublet is increased. 

If the dipole is considered within the framework of the 
discrete element grid where dipole strengths and downwash 
velocities are held constant, or averaged, over each control 
box, a cross-section of the dipole flow pattern would appear 
as in Figure 6b, where the infinite upwash at the doublet 
has been replaced by a finite value w . To determine a 
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value for w q which adequately represents the singularity 

strength within the framework of the discrete element 

approximation, the law of continuity was employed. The 

total fluid outflow in the discrete representation is w q A o , 

where A q is the area of the doublet control box. Modeling 

the entire x-y plane with control boxes without regard to 

wing/wake geometry, the total amount of fluid passing back 

00 

through the x-y plane is Z w A , where w is the average 

n=l n n n 

velocity over the control box with area A n « Therefore, 
continuity requires 



w^A 
o o 



n=l 



w A 
n n 



0 



and if the A ' s are chosen so that A = A , the upwash 
n o n 3 ^ 

velocity is determined by 



w 



o 



Z 

n=l 



w 



n 



( 4 . 9 ) 



In actual practice, the summation of the downwash velocities 
through the x-y plane is necessary only to some finite radius 
from the doublet, since the inverse proportionality of the 
kernel function with distance from the doublet causes the 
value of w Q to converge within a reasonable summation. 

5 . Section and Wing Coefficients 

To obtain section lift and moment coefficients, the 
chordwise pressure distribution at each spanwise station 
was integrated in the manner employed in thin airfoil theory. 
First a coordinate transformation of the chordwise variable 



is made as follows 



u 



CO 




UO 

CT.fc^ 



x /c = ^(l-cos ©) 



The pressure coefficient is then expressed in terms of a 
Fourier expansion of the new variable 

00 

c (0) = c cot % + £ C sin n0 (4.9) 

p p o n=l p n 



where the first term accounts for the leading edge singu- 
larity, while the trailing edge slope singularity is 
acknowledged by the infinite series. The section lift is 
obtained from 
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(4.10) 



Substituting the Fourier expansion for 

coordinate transformation, equation 
c 

p IT 00 

c p = — — / (l+cos0)d0 + E 

2 0 n=l 



0^ and performing the 
(4.10) becomes 

Cp n * 

— — / sin n0 sin0d0 



0 



(4.11) 



Due to the orthogonality of the sine function, this integra- 
tion yields 
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(4.12) 



The section moment about the leading edge is defined by 
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(4.13) 
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Performing the same transformation and integration as for 
the section lift coefficient, the following results 
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(4.14) 



Transferring this moment coefficient to the mid-chord point 
requires the further calculation 




The integration of section lift and moment values to wing 
coefficients is performed in a similar manner. Here the 
coordinate transformation is made to the spanwise variable, 

4 

such that 

y/b = j (1-cos 0) 

The product of the section chord and the section lift 
coefficient is then expressed as the Fourier series. 

00 

0 

cc^ = c. cos — + Z c. sin n0 (4.16) 

o ^ n=l n 



Performing the integration 
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(4.17) 



of the Fourier expansion of the section lift in the trans- 
formed coordinate system, the following results 
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(4.18) 
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The same relationship is obtained for the wing moment 
coefficient when the sectional moment is expressed as in 
the Fourier expansion of equation (4.16). For swept wings, 
the sectional moments were transfered to an axis extending 
from the half-root-chord point, prior to developing the 
Fourier series. 

The infinite series of equations (4.9) and (4.16) are 
of course terminated at the number of chordwise and span- 
wise points respectively of the wing discrete element grid. 
This places a minimum limit on the number of chordwise 
and spanwise control points which are required to achieve 
valid sectional and wing coefficient values, a factor which 
will be discussed in Section VI. 

B. WING/BODY ANALYSIS 

1 . Basic Formulation 

Inclusion of a finite radius body in the nonsteady 
lifting surface analysis is equivalent to adding one more 
boundary condition to the problem formulation: the require- 

ment for no flow through the body surface. For the analysis 
pursued here, the body will be stationary with regard to 
the perturbation motion, and will be idealized as an 
infinitely long cylindrical surface with axis coincident 
with the undisturbed flow. The cantilevered, midmounted 
wing is harmonically oscillating within the limitations of 
small disturbance theory previously developed. 

In steady flow analysis, the body effects have been 
traditionally handled by singularities, matching the wing 
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singularity distribution, at image points within the body in 
the wing plane [39]- This method will satisfy the nonsteady 
boundary conditions only in a quasi-steady sense, because 
of phase differences between disturbances at the body sur- 
face caused by a wing singularity and its corresponding 
image. Therefore, to model the body in the nonsteady 
problem, a system of singularities are placed on the body 
surface establishing a grid similar to that for the wing. 
These curvilinear panels each have an harmonically oscillat- 
ing singularity at its center, coincident with a normalwash 
control point. Thus, in effect the wing grid is merely 
extended over the body surface as shown in Figure 7* 

In order to handle the more complex geometry of the 
wing/body problem, a coordinate transformation from the 

4 

rectangular (x,y,z) system of the wing analysis is made to 
a cylindrical (r,9,x) system, where the undisturbed flow 
direction (x axis) is common. Thus the transformation is 
defined by 



Figure 8 shows the wing/body configuration with the boundary 
conditions in the cylindrical coordinate system. The down- 
wash velocity over the wing is now related to the velocity 
potential distribution by 



x = x 



y = r cos 6 



(4.19) 



z = r sin 0 




0=0 or it 



(4.20) 
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and the no-flow condition for the body surface Is stated as 
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(4.21) 



where r^ is the body radius. As in the wing analysis, the 
wing downwash velocity Uq is determined from the type of 
surface motion prescribed for the problem. 

2 . Governing Matrix Equation 

The governing matrix equation for the wing/body 
problem takes on a much more formidable appearance than the 
basic equation (4.2) for the wing along 
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where {u^} = wing downwash vector 

{u r > = body normalwash vector 

e wing singularity potential vector 

{<f>g} = body singularity potential vector 

The kernel function matrices are defined as follows 



= downwash on wing due to wing singularities 
K^g = downwash on wing due to body singularities 
Kgw = normalwash on body due to wing singularities 
Kg = normalwash on body due to body singularities 

Direct solution of equation (4.22) for the velocity potential 
distribution would be. limited due to the computer storage 
requirements of the large complex kernel function matrix . 
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However this approach is not feasible at any rate, since 
the kernel function matrix is ill conditioned in this form 
and does not lend itself to efficient inversion. 

To attack the problem, equation (4.22) is separated into 
two matrix equations 

*- u e* = [K w H<f> w ) + 

{u^} = [ K gyHby} + [KgHbg) 

The second of these can be solved for the body potential 
distribution 

(* B > - [Kg]' 1 tu r -K BW 4. w ] C-23) 

which is then substituted into the first equation. 

{u 0 } = [K w ] {b w > + tKyg] [Kg] (u r > - [K W g] [Kg] [Kgy] (by} 
This may be rewritten in the following form 

N-KwbV 1 u r } = tWB 1 K BW ]{ V (4 * 24) 

The left hand matrix is a modified wing downwash vector 
incorporating body effects, as does the modified kernel 
function matrix on the right hand side of equation (4.24). 

In the analysis considered here the body is steady at 
zero angle of attack to the main stream, while the wing is 
undergoing nonsteady motion, so that 

{u r > = {0} 

Therefore, equation (4.24) becomes 
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(4.25) 



{u e } 



= [ K w- 



Wb 1 



k bw] 



<v 



which is in the same form as equation (4.2), but in which 
the modified kernel function matrix incorporates the body 
boundary condition. 

3 . Symmetry 

The individual kernel function matrices of equation 
(4.25) incorporating body effects can be quite large when 
control points are placed around the entire circumference 
of the body, severely limiting the utility of this method. 
However these matrices may be reduced appreciably in size 
through use of symmetry. Considering the representation of 
Figure 9a, it can be seen that for wing motion symmetric 
about the x-z plane 




} " { *W 



1 



}; 




} 

2 




}; 




} 

3 



U R ) (4.26) 
B 4 



where (<f>^ } stands for the doublet strength distribution on 
i 

the respective wing, and {<f>g } stands for the singularity 

l 

strength distribution on the body in the respective quadrant. 
It can also be seen that due to the antisymmetric nature of 
the wing doublet flow with respect to the x-y plane 



{(}> } = - {A }; {<f> } = - (4> } (4.27) 

ttq b 3 B 2 



since the purpose of the body singularities is to counter 
the flow of the wing doublets through the body surface. 

Considering first the body kernel function matrix (Kg), 
which is defined by the relation 
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{u r } = [K b ]( 4> b } 



or 



(u } = [K J(<J> R } + [K R m R } + [K r ] {<t) R } + [K r ]{<|) } 

r B 1 B^. b 2 B 2 B 3 B 3 B^ 

Symmetry allows this formulation to be reduced to 



( u r > [K b 1 + K B 2 K B 3 K B 1 | I ^B 1 ^ 



(it. 28) 



In the case of the wing-body kernel function matrix (Kyg)j 
this same consideration produces 



{u e } = 



[K 



WB. 



+ K 



WB, 



- K t . id - K t ., d ]{<t> D ) 



WB 



WB, 



B-, 



(4.29) 



Control points need only be placed, therefore, over the 
first quadrant of the body surface, reducing the size of 
the body effect kernel function matrices of equation (4.25) 
by one-fourth. The surface over which the flow conditions 
are specifically satisfied are indicated by the solid line 
of Figure 9b, which then satisfies the boundary conditions 
on the whole wing/body surface through symmetry. The body- 
wing kernel function matrix (Kg^) follows the same develop- 
ment as in Section IV. A. 3 for the wing alone and has the 
same form as given by equation (4.8). 

It should also be noted that for antisymmetric wing 
motion the kernel functions of equations (4.28) and (4.29) 
would have the form 



[K] = [K 1 - K 2 + K 3 - K h ] 
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Since any general harmonic motion of the wing can be 
expressed as a combination of symmetric and antisymmetric 
motion, this method, as in the wing alone case, need not 
be restricted to the symmetric motion case considered in 
detail here. 

4. Effective Body Length 

The body considered in this analysis is idealized 
as an infinitely long cylinder with axis aligned with the 
undisturbed free stream flow. This is not an unrealistic 
limitation, since in most practical cases, the center part 
of a fuselage is nearly cylindrical and the fineness ratio 
of the fuselage is large enough so that the flow at the 
center part is nearly the same as for an infinitely long 
cylindrical body. Experiments in steady flow have proved 
that beyond an effective body length, the lift distribution 
of the wing/body combination is independent of the body 
length [ 4 0 ] . Steady flow analyses, such as Woodward's [ 4 1 ] , 
employ a "wing-body interference region," where the body 
no-flow boundary condition is explicitly applied a finite 
distance upstream and downstream from the wing root section. 

A similar effective body length is modeled in this 
approach, where the body singularity grid extends from a 
finite distance upstream of the wing root to a point down- 
stream of the effective wake. The actual body length which 
must be modeled in order to include all interference effects 
is discussed in Section VI. Since the kernel function . 
velocities decrease approximately as the distance from the 
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singularity cubed, it is not unreasonable to assume that 
interference effects are concentrated close to the wing/wake 
area. 

5. Body Singularities 

The acoustic singularity used to model the wing was 
a potential doublet, due to the requirement for a pressure 
differential across the wing surface. No such requirement 
exists for the body singularities, so that an harmonically 
oscillating source could be employed as well as the dipole. 
However, in the development of the wing computer program, 
it became apparent that the solution to the matrix equation 
(4.4) was very sensitive to the value placed on the kernel 
function upwash singularity. Approximate methods used to 

obtain an average value of the upwash did not produce 

* 

acceptable loading distributions over the wing. Only when 
this upwash was determined numerically, as indicated in 
Section IV. A. 4 were good results obtained. Unfortunately 
a source singularity, while having a somewhat simpler form 
of kernel function, can not be analyzed by this same type 
of numerical development. 

The body was therefore modeled by a network of harmon- 
ically oscillating potential doublets, as on the wing, with 
axes oriented in the radial direction. Each doublet is at 
the center of its control panel, coincident with a normal- 
wash control point. Determination of the normalwash (radial 

velocity u ) at the doublet location is obtained, as with 
r 0 

the wing singularity, from a consideration of continuity. 
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The total fluid outflow from the doublet is u^ A q , where 

A q is the area of the doublet control panel on the body. 

The amount of fluid passing back through the body at a 

control panel away from the doublet is (u ) A where n 

r nm nm 

indicates the axial location and m the circumferential loca- 
tion of this control panel. Continuity is then expressed 
for this discrete element representation as 



oo m 

u A + E E 
r o , , 

o n=l m=l 



( u ) A, 



r nm nm 



0 



where M equals the number of panels located around the 

circumference of the body. If the A ' s are chosen so that 

nm 

A nm = A q the value of the doublet singularity is given by 

oo m 

u = - E E (u ) < (4.30) 

o n=l m=l nm 



This continuity equation states, in effect, that all the 
fluid which leaves the doublet in the radial direction must 
pass back through the body surface prior to returning to 
the doublet. As in the wing singularity case, the summation 
of normalwash in the axial direction can be terminated after 
a reasonable distance in both the upstream and downstream 
directions . 

6 . Kernel Functions 

The elements of each of the kernel function matrices 
defined in Section IV. B. 2 represent the normalwash on either 
wing or body surface due to a harmonically oscillating 
doublet of unit strength which satisfies the governing 
linearized differential equation (2.12). Doublets on the 
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wing are of course oriented in the 0 (or z) direction, 
while those on the body have their axes in the radial 
direction . 

The formulation of these kernel functions is developed 
in Appendix A. They are summarized below with control 
point at (r,0,x) and doublet located at (r ,0 o ,x o ). The 
body radius is r^. 

(i) - Downwash on wing due to unit wing/wake singu- 

larity . 

2 

u fi = — a (1 + i — R)exp {i — ? [m (x-x )-r| } (4.31) 

6 4irR J a 6 a 3 L ° J 

00 00 



where 



•R = -x/(x-x o ) 2 + 6 2 (r 2 +r^ + 2rr Q ) , 



(ii) - Normalwash on body due to unit wing/wake 

singularity . 



“r- 6 + 1 7T5 R) 



4ttR' 



a B 

oo 



+ { Te? ) t 6 " ^ exp u 772 (1, - 32:i 






where 



R = \/(x-x o ) 2 + 6 2 (r 2 +r 2 + 2r Q r B cos 0) 



3R _ B‘ 



3r r B - o 



(r c + r cos 0) 
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In both the above formulations, the minus sign represents a 
singularity on the 0=0 wing/wake, and the plus sign a 
singularity on the 0= tt virtual wing/wake. 

(iii) Kg - Normalwash on body due to unit body singu- 
larity . 



u = -=& 



4ttR : 



cos(0-0„) + ^ r B (1 - cos<0-0 Q >) 



O' 



1 + i 



exp \ i 



■]- 

l 1 rts [ M » (5C - x o > - R J) 



a (3 

oo H 



(— ^) Rr (1- cos<0-0» 15 
a B 2 B o 3r 



(4.33) 



a B‘ 

oo H 



where 



R = + 2B 2 rg (l-cos<0-e o » 

I? = T r B (1 -°° s < e - 0 o >:> 



(iv) Kyg - Downwash on wing due to unit body singu- 
larity . 

,2 



U 0 = 



-B 



4ttR : 



„ i 3 t „ , 3R 

r sin 6 0 + r (r B - r cos 0 Q ) 






i — r| - ( - ? ) R(r R - r cos 9«) 9R 

aj 2 J a ro B 2 B 



V 30 



where 



exp (i [m(x-x ) - r]) 

l L ° JJ 



(4.34) 



R 



= \/(x-x o ) 2 + B 2 (r 2 +r B - 2rr B cos 0 Q ) 



3R _ B 

30 - " TT rr B sin 0 o 
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V. DESCRIPTION OF COMPUTER PROGRAMS 



A. GENERAL DESCRIPTION 

Listings of the two computer programs developed and used 
in this investigation are reproduced in Appendices D and E. 
The programs are written in standard FORTRAN IV language. 
Calculations were performed on the I.B.M. 360/67 computer 
at the W. R. Church Computer Center of the Naval Postgraduate 
School. Object codes were obtained with the I.B.M. G-level 
compiler. 

Both the wing and the wing/body programs are arranged in 
the same general format. The MAIN program reads the input 
data and establishes the wing or wing/body control grid. 

4 

One or more subroutines are called which establish the 
kernel function matrix elements. MAIN calculates the wing 
downwash velocity vector and calls subroutine COMAT to 
solve the matrix downwash equation (4.2) for the wing velo- 
city potential distribution. MAIN finally calls subroutine 
PRES which calculates the perturbation pressure distribution 
over the wing by applying finite difference approximations 
to equation (4.5). PRES then calls subroutines SECLM and 
WINGLM to integrate this distribution to obtain sectional 
and wing forces and moments . These results are printed by 
PRES and control is returned to MAIN for the reading of 
data cards for the next problem. 
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B. WING PROGRAM 



The wing program computes potential and pressure 
loadings due to harmonic motion on general planar wing 
configurations from rectangular to arbitrary sweep angles 
for both the leading and trailing edges. The program is 
restricted, however, to constant sweep angles, and to 
a finite tip chord with minimum taper ratio of about one- 
fifth. Thus delta wing configurations are excluded. This 
restriction is caused by the method used to integrate the 
chordwise pressure distribution to obtain sectional lift 
and pitching moment, discussed in Section IV. A. 5* The 
program provides for a maximum of 100 control points on 
the wing and ten control points in either the chordwise 
or spanwise directions. Storage requirements are the 

4 

equivalent of 18,000 single precision complex words, or 
144,000 bytes on the I.B.M. 360, with a maximum run time 
of approximately twelve minutes for one problem. Provision 
is made for the running of successive problems for as many 
sets of input data cards as are provided. 

Figure 10 is a diagram of the wing program, while sub- 
routine flow diagrams are presented in Appendix C. Sub- 
routine DINCO forms the kernel function matrix from 
equation (4.8), with the matrix elements defined by equation 
(4.3). The COMAT subroutine solves the linear complex 
matrix equation (4.2) by the Gauss-Jordan method with total 
pivoting. This subroutine was written for systems of real 
equations by Mrs. Sharon Good, David Taylor Model Basin, 
as published in Ref. 42, and modified to include the complex 
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capability by Mr. Hellmut Golde, Department of Electrical 
Engineering, University of Washington. COMAT was further 
modified by the author for particular application to the 
wing and wing/body programs. 

Input instructions are presented in Appendix B. Any 
number of spanwise and chordwise control points, up to the 
maximum of ten each, may be specified for a semi-span wing. 
The wing is modeled by an equal number of chordwise control 
points at each spanwise station. The program is limited 
by theory to the subsonic flow regime, within which any 
mode, amplitude, or frequency of harmonic wing motion may 
be specified, including the steady case. Results printed 
by the program for each spanwise station are: 

(i) Potential doublet strength distribution 

4 

(ii) Pressure coefficient distribution 

(iii) Section lift and pitching moment coefficients 

In addition, wing lift and pitching moment coefficients are 
presented. 

C. WING/BODY PROGRAM 

The wing/body program incorporates the body surface 
boundary condition into the harmonic motion analysis of a 
rectangular wing planform. As in the wing program, a maxi- 
mum of 100 control points may be used to model the wing, 
while up to 130 points are allowed in the body control panel 
network. Storage requirements are the equivalent of *13 >500 
single precision complex words, or 3^8,000 bytes on the 
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I.B.M. 360. A maximum run time of about 20 minutes for 
the wing/body problem is required, while a wing only solu- 
tion (body radius equals zero) is obtained in less than 
seven minutes. As in the wing program, successive problems 
for different wing, body, and flow geometries may be run. 

Figure 11 diagrams the wing/body program, with indivi- 
dual subroutine flow charts presented in Appendix C. The 
influence coefficient matrices defined in Section IV. B. 2 
are formed in the following subroutines: 

(i) DINCO - Dy from equation (4.31) 

(ii) UTHETA - from equation (4.34) 

(iii) URAD1 - Kgy from equation (4.32) 

(iv) URAD2 - Kg from equation (4.33) 

The matrix K is inverted by subroutine COMAT and the modi- 
B 

fied kernel function matrix of equation (4.25) is formed 
in MAIN. COMAT is again called to solve the resulting matrix 
equation and PRES performs the same function with the same 
data printouts as in the wing program. Subroutine DINCO, 
while performing the same function as in the wing program, 
is here in a more simplified form because of the restricted 
wing geometry of the wing/body program. This accounts for 
the reduced running time of wing only problems in this 
program as compared to the general wing program. 

The wing/body program uses the same wing grid and flow 
geometry parameters as the wing program. Any number of 
body control points may be specified up to the maximum of 
130. An effective body 'length of one wing chord upstream 
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from the wing leading edge to one wing chord downstream 
from the termination of the effective wake is modeled by 
the program. Input data instructions are presented in 
Appendix B. 



VI. RESULTS AND DISCUSSION 



A. WING ANALYSIS 

1 . Comparison with Lifting - Surface Theory 

In order to investigate the validity of the discrete 
potential element approach to lifting surface theory, the 
wing program was run for a wide variety of wing/flow geometries 
and types of oscillatory motion. Examples were picked which 
could be checked against previous work reported in the 
literature, representing a variety of lifting surface theory 
approaches, as well as the relatively limited amount of 
experimental results which have been obtained in the unsteady 
field. 

4 

The figures presented in this section are a representative 
summary of this investigation. The subsonic flow regime was 
covered from the incompressible M co =0 to the high subsonic 
M oo =0.9, while the range of frequencies varied from the steady 
case to the relatively high reduced frequency of k= 1.2. 

The wing planforms considered had a range of aspect ratios 
from one to six, sweep angles from zero to 45 degrees, and 
taper ratios from one to one-half. The types of symmetric 
harmonic motion considered were: 

(i) Pitching - rigid body oscillations of the wing 
about a spanwise axis perpendicular to the flow direction; 

(ii) Bending - oscillations of the wing as a beam canti- 
levered at the root chord in the first natural mode of bending. 
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(iii) Plunging - rigid body oscillations of the wing in 
the z direction at a mean angle of attack of zero; 

(iv) Flapping - rigid body roll oscillation of the wing 
simply supported from a chordwise axis inboard of the root. 
Steady data were obtained with the wing analyzed as a rigid 
body at a fixed angle of attack. The results have been pre- 
sented as chordwise pressure distributions, spanwise lift 
and pitching moment distributions, and wing coefficients 
plotted with respect to reduced frequency. 

In each of the figures, the coefficient amplitudes are 
normalized with respect to the motion as follows: 

(i) Pitching - angle of attack amplitude; 

(ii) Bending - wing tip angle of attack amplitude; 

(iii) Plunging - vertical motion amplitude; 

4 

(iv) Flapping - amplitude of flapping angle. 

These represent the conventions used in the applicable 
references , however the phase angle relationship were con- 
verted to the coordinate convention used in this paper where 
differences occurred. 

Correlation of the wing program results with existing 
lifting surface theory was in general quite good. Devia- 
tions appeared to come from the different methods employed 
to handle the wing edge singularities in the pressure 
distribution. The lifting surface approaches assume the 
form of these singularities a priori in the choice of their 
loading functions, as discussed in Section III.B. The 
discrete element approach, on the other hand, assumes no 



loading profile , but obtains proper potential and pressure 
distributions through inclusion of the boundary conditions 
in the problem formulation. As discussed in Ref. 1, much 
effort has been devoted in lifting surface theory to 
developing the most numerically efficient and accurate 
loading functions, as well as to improve the numerical 
methods of handling the singular kernel functions. 

The first four figures compare wing program results with 
one of the most recent lifting surface theories. This 
advanced kernel function method, based on the acceleration 
potential, was presented by Laschka in 1963 [5] and further 
developed by Laschka and Schmid for interfering planar 
surfaces in a 1967 paper [43]. Figure 12 compares chordwise 
pressure distributions on a swept, tapered wing undergoing 

4 

pitching oscillations in low subsonic flow with Laschka' s 
results [46]. Good correlation is evident at each spanwise 
station, with the only variance being a small difference in 
the shape of the pressure distribution near the leading 
edge, as previously discussed. 

Figures 13 and 14 present spanwise loading on a 45 degree 
swept constant chord wing in both the pitching and plunging 
modes. Results are compared with Lashka's work [5] for both 
low subsonic (M =0) and M =0.8 flows. Correlation is 
again good. The discrete potential element approach appears 
to overestimate the amplitude of the lift and pitching 
moment slightly, especially near the wing tip. In fact this 
latter trend appears in almost all the results of the present 
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method. The explanation again would appear to lie in the 
method of handling the wing tip pressure slope singularity, 
which is included implicitly in the boundary conditions. 
Further development of the discrete potential element 
approach would definitely have to include investigation of 
the adequate recognition of this wing tip condition in the 
problem formulation. 

In Ref. 37 Laschka compares his results, for the wing 
configurations of Figures 13 and 14, with the work of Pao 
Tan Hsu [44] . It should be noted that the wing program 
results agree more closely with Laschka' s data, while fall- 
ing between Laschka' s and Hsu's results except for amplitudes 
at the wing tip. The theory here is further compared with 
experimental results obtained by Laidlaw [45] with fair 
correlation . 

Figure 15 compares wing program spanwise loading with 
Laschka' s results from Ref. 37 for a steady swept, tapered 
wing in low subsonic flow. The lift and pitching moment 
coefficients are presented on an expanded scale, with the 
same type of comparisons already noted. 

The results presented in Figures 16 through 18 are note- 
worthy for the Mach number range considered (M oo =0.24 to 0.9) 
and for the consistent experimental data presented from the 
work of Lessing, Troutman, and Menees [47]. Spanwise and 
chordwise loadings are plotted for an aspect ratio three 
rectangular wing in the bending mode. Theoretical comparison 
is made to the lifting surface results developed by Lessing 
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et al in I960 from the original kernel function methods of 
Refs. 25 and 48. 

At a Mach number of 0.24, the chordwise pressure loadings 
obtained from the two theories are almost identical. The 
spanwise lift and moment curves show a somewhat heavier 
loading concentration near the wing tip for the discrete 
potential element method, as in the preceding figures. This 
trend also appears at Mach numbers of 0.7 and 0.9* The 
M^ = 0.7 chordwise pressure distributions show a difference 
in the representations of the leading edge pressure distri- 
butions outboard of n = 0.5. These deviations can probably 
again be laid to the different methods employed to handle 
the wing pressure singularities. However, the wing program 
results agree more closely to the experimental data of Ref. 
47, than does the kernel function approach. No explanation 
can be found for the deviation of the pressure phase angles 
in the wing trailing edge area for both the 0.24 and 0.7 
Mach number cases. It should be noted, however, that the 
experimental Doints again correlate more closely with the 
wing program results in this area. In fact, the close cor- 
relation of these experimental data with the results of the 
present method is quite gratifying. 

It Is interesting to observe that a shock wave may have 
started to form on the wing in the 0.7 Mach number flow down- 
stream of the 60 % chord line, as indicated by the drastic jump 
in the pressure phase angle measurements at this location in 
Figure 17b. The theoretical analyses would not, of course, 
reflect the existence of the shock wave. 
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The closeness of the two theoretical approaches at 
M^ = 0.9 is also noteworthy. As discussed in Section II. C, 
the linearization of the basic problem for perturbation 
analysis would appear somewhat questionable this close to 
the speed of sound. However, the consistency of the results 
would seem to indicate that the linearized subsonic theory 
is still valid at this high subsonic Mach number. 

Figure 19 analyzes the pitching motion of an aspect 
ratio two rectangular wing. Chordwise and spanwise loadings 
are compared with the experimental and theoretical results 
of Laidlaw [45]. The latter represents a numerical treat- 
ment of the rectangular wing aspect ratio theory of Reissner 
[49, 50] developed in 1947- Considering the fact that 
Reissner' s approach basically involves applying correction 

4 

factors to two-dimensional results in order to account for 
finite aspect ratio, this theory agrees quite well with the 
discrete potential element approach. The experimental data, 
although not as consistant as that of Lessing [47], agrees 
reasonably well with the Wing Program results. Figure 20 
presents spanwise loading for the same wing in plunging 
motion. Correlation of pressure distributions between 
theories and experimental data for this mode of nonsteady 
motion is the same as for the pitching case. 

The frequency response of a rectangular wing in the flap- 
ping mode is presented in Figure 21. Comparison is made to 
experimental and theoretical data by Woolston, Clevenson, 
and Leadbetter [51] , who employed a basic kernel function 
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approach. Correlation of the wing program results with the 
experimental data is reasonably good, although there is some 
difference in the phase angle values. However, the varia- 
tions of both lift and pitching moment with reduced fre- 
quency compares very well. Woolston's theoretical points 
coincide essentially with the wing program curves except 
for the pitching moment phase angle. No explanation can be 
found for this difference. Laschka [5 ] > in comparing his 
results with those of Woolston, produced virtually the same 
curves as the discrete potential element program. 



% 
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FIGURE \4b 
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FIGURE I7 to 
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2 . Convergence 



Figures 22 through 24 present spanwise loading for 
three different wing planforms as a function of the singu- 
larity grid density in order to show convergence of the wing 
program. Three types of motion are considered, with wing 
grids numbering from 25 to 100 points. Figure 25 shows 
convergence of the wing program for the steady case on an 
aspect ratio six rectangular wing. This wing planform is 
used in the basic reference for the wing/body analysis in 
steady flow. The ordinate of these graphs has been greatly 
expanded to permit a good evaluation of convergence. 

All of the examples are seen to converge at about the 
same rate. At 50 control points on the wing the maximum 

deviation is less than four per cent. At 60 or more points 

« 

the results are virtually coincident. The shape of the 
control box would appear not to be a factor, since varying 
the grids on the individual wings of Figures 22 through 25 
also varies the shape of the individual boxes. However, it 
would seem prudent to maintain the control boxes at reason- 
able aspect ratios (less than about ten) to assure valid 
results. Another restriction on the wing grid size is the 
requirement to integrate the pressure loading to obtain 
section and wing coefficients. A minimum of six points in 
either the chordwise or spanwise directions was found neces- 
sary to adequately define the loading distribution for the 
integration subroutines. However, the basic requirement for 
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6 0 or more points in the wing grid makes this latter 
restriction a factor only for higher aspect ratio wings. 

A second type of restriction on the wing grid size is 
caused by the accuracy of the inversion method of subroutine 
COMAT. Solution of the downwash matrix equation is accom- 
pli shed by the Gauss-Jordan method employing total pivoting. 
With the relative magnitudes of the kernel function matrix 
elements occurring in this type of analysis, the inversion 
process starts to lose accuracy with a system of equations 
or order greater than about 110. Computational accuracy 
was verified by substituting the velocity potential vector 
solution into the downwash matrix equation, and obtaining 
a downwash vector to compare with the original input. Checks 
of the solution for grid densities up to 150 points were 
made, for various wing/flow conditions, with consistent loss 
of accuracy. Thus, the program has been limited to wing 
grids of 100 control points or less to ensure accuracy. 

This also limits the aspect ratio of the wings that may be 
analyzed; however, for the purposes of this report, aspect 
ratios of six or less were well within the capability of 
the program. If added capability were required, a more 
sophisticated inversion routine would have to be developed. 
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FIGURE Z2a 
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FIGURE Z3b 
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FIGURE 24 a 

CONVERGENCE - FLAPPING 



k= 04 AR = 3 




104 



FIGURE ZAb 
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3. Wake Effect 



Figures 2 6 and 27 summarize the effect of the finite 
wake length considered in the singularity grid network as 
discussed in Section IV. A. 2. Seven configurations of wing/ 
flow geometries were analyzed, covering the range of condi- 
tions considered in this paper. Wake effects were included 
to a maximum distance of ten root chord lengths downstream 
from the wing trailing edge. In all cases, wing loading 
had converged to within one per cent of amplitude and one- 
half degree of phase angle by four root chord lengths of 
effective wake. This correlates with Haviland's work [8], 
in which he reported that the effect of the wake on the wing, 
for a rectangular planform in the steady case, was determined 
to within one per cent in five chord lengths . 

4 

Figure 26 presents a specific case of wing coefficient 
variation with effective wake length. The maximum deviation 
of the wing coefficients, with respect to an effective wake 
length of four root chords, is shown in Figure 27. From 
the latter figure, it can be seen that the near wake (within 
one root chord length of the trailing edge) has a very great 
effect on the wing loading. The wake influence then dies 
out quite rapidly as the effective termination point is 
moved downstream. Therefore, the finite wake, taken into 
account in determining the aerodynamic loading on the wing, 
need only extend a realtively short distance downstream in 
order to obtain accurate results. 

In the wing analysis of this paper, an effective wake of 
four root chord lengths was used. The wing program will. 
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however, incorporate any effective wake length desired by 
varying an input parameter (Appendix B). With a wing grid 
of 100 control points, four root chord lengths of effective 
wake requires a wake grid which varies from 400 singularities 
for a rectangular wing, to approximately 1100 singularities 
for a wing with a taper ratio of one-fourth. The increased 
number of wake singularities are required for a tapered wing 
because of the decreasing chordwise dimension of the control 
boxes towards the wing tip. Thus a greater number of 
singularities are required, at a spanwise station outboard 
of the root, to extend four root chord lengths into the wake, 
than are required for a rectangular wing. The large number 
of wake singularities pose no computational problem, since 
the effect of the wake is incorporated prior to the solution 

4 

of the kernel function matrix equation, as discussed in 
Section IV. A. 2. 
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PER CENT DEVlAT\ON IN AMPLITUDE 
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B. WING/BODY ANALYSIS 



1 . Comparison with Steady Analyses 

A thorough literature search was made in an attempt 
to find previous theoretical or experimental work, in the 
nonsteady lifting surface field, with which to compare and 
validate the wing/body program results. Unfortunately, no 
such effort which explicitly analyzed body interference 
effects on an oscillating wing's pressure distribution could 
be found. Rodden [38] includes some body interference effects 
in his doublet lattice method for analyzing nonplanar con- 
figurations, by extending the lifting surface elements onto 
the body surface near the wing-body intersection. However, 
the actual interference effects are not presented. 

It was therefore necessary to compare the wing/body 
program results with data previously obtained in the steady 
field. The two latest works found were presented at the 
AGARD conference on Aero-dynamic Interference in 1970 by 
Kuchemann [52] and Labrujere [53] • Both analyzed an aspect 
ratio six rectangular wing at a six degree angle of attack 
midmounted on a cylindrical body of radius approximately 
equal to 20% of the wing semi-span. This configuration was 
also theoretically analyzed by Weber [40] and experimentally 
investigated by Korner [54] . Each of the theoretical analyses 
employed a combination of surface source singularities and 
vortex line distributions to construct a steady kernel func- 
tion method solution. - 

Figure 28 compares the wing/body program results, for 
both the wing and wing/body combination, with results 
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obtained from the above investigations. The discrete 
potential element approach spanwise loadings are seen to 
compare almost exactly with the theoretical curves obtained 
by Weber, while the experimental data of Labrujere and 
Korner offer good correlation. The theoretical work of 
Labrujere overestimates the spanwise loading, especially in 
the case of the wing/body combination. The good agreement 
between these theoretical and experimental results and the 
wing/body program results provides a reasonable measure of 
confidence in the discrete potential element analysis of 
this type of nonplanar configuration. 

Figure 29 indicates the convergence of the wing/body 
program as the number of control points on the quarter cir- 
cumference of the body was increased from 72 to 108. This 
graph is a greatly expanded representation of the wing root 
area effects. Outboard of the 25% semi-span point, the 
three grid configurations gave essentially the same results 
shown in Figure 28. In each case, body and wing control 
panels were matched as closely as practical in size and 
shape, so that the body grid becomes in essence an extension 
of the wing network on the body surface. Small variations 
between wing and body panel size did not effect the wing 
pressure distribution; however, when wing and body grids 
were obviously mismatched results became inconsistent. 

For all the data presented in this paper, the body grid 
extended from one chord length upstream of the wing leading 
edge to one chord downstream of the effective wake termina- 
tion point. Extending the body grid somewhat further in 
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either direction produced no appreciable effect on the 
results. However, decreasing the length of the body grid, 
to the point where the wing/wake singularity grid extended 
beyond the body grid, caused the wing pressure distributions 
to become inconsistent. This effect held for both the 
steady and nonsteady cases. Therefore, it would seem that 
the body surface which effects the wing pressure distribu- 
tion, as discussed in Section IV. B. 3, is concentrated in 
this effective length between one wing chord upstream and 
downstream of the wing/wake grid. Further investigation of 
this effect, such as by increasing the effective wake and/ 
or the body grid length with increased number of control 
points, was precluded due to the accuracy limitations of 
the matrix inversion routine discussed in Section VI. A. 2. 

4 

In addition, wing body grids within this effective body 
length were limited in the number of control points which 
could be used for the same reason. It was, therfore, not 
possible to investigate configurations with body radii 
greater than about 20% of the wing semi-span. Future work 
with this method would definitely require an improved com- 
plex matrix inversion procedure to allow greater range in 
the wing/body analysis. 

An indication of the effect of body size, in the steady 
case, is presented in Figure 30, where spanwise loadings 
are plotted, with a greatly expanded ordinate, for the wing 
previously considered -in combination with bodies of three 
different radii. It can be seen that the interference 
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effect is large even for a relatively small body radius of 
five per cent of the wing semi-span. This is reasonable 
since the body is at zero angle of attack and is providing 
no lift carry-through, or reinforcement, for the wing. This 
effect for bodies of small radius is also shown by Kuchemann's 
results [52] . 
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FIGURE 28 
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2 . Extension to Nonsteady Interference 



Investigation of wing/body interference for the non- 
steady case was conducted for two modes of wing motion: 
bending, as previously discussed; and torsional oscillations 
of the wing cantilevered at the root in an assumed first 
natural mode. In each case the coefficients were normalized 
with respect to the wing tip angle of attack. It should 
be noted that, of the types of motion normally considered 
in the nonsteady case, these are the only two valid for the 
wing/body configuration considered in this investigation. 

Figure 31 presents the interference effects for the wing 
analyzed in Figure 16, cantilevered from a body with ten per 
cent semi-span radius. The decrease in sectional lift 
amplitude follows the steady case format, being greatest at 

4 

the root and essentially disappearing as the wing tip is 
approached. The change in phase angle, while small, remains 
relatively constant until midspan and then slowly decreases. 
Figure 32 presents spanwise loading in the root area for 
wing/body configurations with body radii varying from five 
to twenty per cent wing semi-span. Again the trends are 
essentially the same as in the steady case, with a small but 
relatively constant phase angle difference. Wing/body inter- 
ference for the torsional vibration case, with body radius 
equal to 20% of the wing semi-span is shown in Figure 33* 
These interference effects in the torsional loading follow 
the same trends as for the nonsteady bending mode, as well 
as the steady case. 
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The results of the wing/body program show that nonsteady 
interference occurs in the same way and via the same mech- 
anism as in the steady case. Interference is greatest in 
the root area, decreasing towards the wing tip. The dif- 
ference between the steady and nonsteady cases comes from 
the type of wing loading with which the body interfers. For 
a wing at a steady angle of attack, the wing loading is 
greatest in the root area, and, therefore, the interference 
is of relatively large magnitude. In the nonsteady case, 
the wing loading is smaller in the root area (since the wing 
is cantilevered from the body and has no notion at the root) 
and increases towards the tip, as the wing motion amplitude 
increases. Therefore, the interference effects, concentrated 
towards the root area, are of smaller relative magnitude for 
nonsteady motion, as compared to the general steady case. 
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FIGURE 3Zb 
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VII. CONCLUSIONS 



In summary , a discrete potential element approach to 
subsonic numerical lifting surface theory has been developed 
and shown to be practical in predicting the nonsteady load- 
ing on harmonically oscillating wings. This approach was 
then extended to the case of an oscillating wing, canti- 
levered from a steady cylindrical body, to investigate inter- 
ference effects and show the versatility of the basic method. 

Correlation of the wing program results with those of 
nonsteady lifting surface theory is generally quite good 
over the full range of subsonic flow. Deviations, where 
they exist, appear to come from the different methods 
employed to handle the wing edge singularities in the pres- 
sure distribution, which are assumed a priori in the lifting 
surface theory, but in the discrete potential element 
approach are implicit in the boundary conditions. Primary 
deviation appears to rest in the handling of the wing tip 
pressure slope singularity, since sectional lift and pitching 
moment are generally overestimated in the wing tip area. 
Future work with the discrete potential element method should 
include investigation of a more adequate acknowledgement of 
this boundary condition in the problem formulation. 

A minimum of about sixty control points is required on 
the semi-span wing surface to achieve convergence of the 
wing program results. In addition, at least six control 



126 



points are required, in both the chordwise and spanwise 
directions, to achieve accurate integration of the pressure 
distribution into sectional and wing coefficients. The 
wings, which can be analyzed efficiently by this method, are 
effectively limited to medium to small aspect ratios, because 
of the requirement for large numbers of control points, and 
consequently long computer run times, for high aspect ratio 
wings. The wing is further restricted to constant leading 
and trailing edge sweeps, and to a finite tip chord. Future 
development of the program capability should be pointed at 
removing these restrictions. 

The effect of the wake on the wing loading is concen- 
trated in the area just downstream of the wing. In all cases 
investigated, wing loading converged to within one per cent 
in an effective wake length of four root chords downstream 
of the trailing edge, allowing the wake singularity grid to 
be terminated at this finite distance from the wing. In 
addition, the wake effect is included in the wing kernel 
function matrix prior to the inversion process, through 
application of the boundary conditions in the wake. In 
this way, an historical drawback to the velocity potential 
formulation, the requirement to integrate the lifting 
surface downwash equation over both wing and wake , is 
removed. 

Wing/body program results agree very well with both 
theoretical and experimental results for the steady case. 
Results obtained for nonsteady wing motion appear consistent 
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and give a good representation of interference effects. 
Convergence of the wing/body program is achieved with a body 
control panel grid of essentially the same format as the 
wing grid. 

The effective body length which causes interference in 
the wing pressure loading extends approximately from one 
chord length upstream of the wing leading edge to one chord 
length downstream of the effective wake termination point. 
Increasing the body length modeled did not effect the 
pressure distribution, while decreasing this length so 
that the body grid did not extend beyond the wing/wake 
singularity grid in either the upstream or downstream direc- 
tions caused inconsistent results to occur. This analysis 
was somewhat limited due to numerical restrictions on the 
allowable size of the body grid. 

In the numerical analysis, the number of control points 
on the wing and the body was limited due to loss of compu- 
tational accuracy in the matrix inversion routine of sub- 
routine COMAT. This restricted the scope of the wing/body 
interference investigation both as to the length and radius 
of body which could be considered. Accurate solution of 
the wing downwash equation and inversion of the body influence 
coefficient matrix, is limited to coefficient matrices of 
order less than 110. Future wing/body analyses by this 
method would definitely require a more sophisticated inver- 
sion routine capable of handling much larger body grids. 
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Within the limitations noted, the wing/body program 
results show that interference effects are significant, and 
follow the same format, in both the steady and nonsteady 
cases. The decrease in wing pressure loading amplitude, 
caused by the body's presence, is greatest in the root area, 
decreasing towards the tip. Differences in phase angle, 
while small, exist over the entire wing. Therefore, these 
interference effects are of importance to three-dimensional 
analyses of wing/body configurations. 

Future development of the discrete potential element 
approach should include the analysis of oscillating wings 
with control surfaces. This approach appears to provide a 
direct means of including the control surface boundary con- 
ditions, without the problems associated with properly 
defining singularity effects in the continuous loading 
formulation. 
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APPENDIX A 



KERNEL FUNCTION DEVELOPMENT 

For the purposes of nonsteady lifting surface analysis, 
a kernel function is defined as the normalwash on a surface 
due to a unit oscillatory acoustic singularity which 
satisfies the linearized convective wave equation (2.9). 

The singularities may be elementry radiators of zero order 
representing simple point sources, commonly used in super- 
sonic analyses, or first order radiators, namely dipoles or 
doublets, whose axes are normal to the surface, used in 
subsonic analyses as in this paper. Development of the 
general form of the doublet kernel function was first 
accomplished by Kussner in 19^0 [2]. The development in 
this appendix follows Kussner 's method and specializes the 
results to the forms of the kernel function employed in 
the discrete potential element approach for both the wing 
and body surfaces. 

In summary, the solution to the linearized wave equation 
for a stationary acoustic singularity is first extended to 
the case of a moving singularity in a space fixed coordinate 
system through application of a Lorentz coordinate transfor- 
mation invariant with respect to the speed of sound. This 
solution is then transferred to the moving wing (or body) 
fixed coordinate system through a Galilean transformation. 
The results are subsequently reduced for the harmonic case 



to the forms of the kernel functions necessary in the analy- 
sis discussed in Section IV. 

1. Stationary Singularity 

The perturbation analysis of Section II led to the 
linearized convective wave equation (2.9), which for the 
stationary case (U^ =0) has the familiar form 



Since the pressure and velocity potential are linearly 
related by 

P = " P ro ^ ( A - 2 

at 

the purturbation pressure also satisfies the wave equation 



Considering a stationary source at the origin of the coor- 
dinate system, equation (A. 3) becomes 




(A. 1) 




(A. 3) 
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(A. 4) 



where 




This equation has the well known solution 



p = ~ F (t-r/aj + £ G (t+r/aj 



(A. 5 ) 
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where F and G are arbitrary functions. From the boundary 
conditions of the problem, only the solution representing 
outgoing waves (F) is admitted. 

Representing the oscillating source as a sphere of radius 
r fe expanding and contracting with a rate of flow away from 
the surface defined by 



The momentum equation (2.2) becomes in linearized form 



S(t) = ilTTr^ u r (t) 



1 3F _ F_ 

r dr 2 
r 




8u 

r 



or 



1 3F _ _F 

r 3r 2 
r 




b 



(A. 6) 



If r b is very small, -pj- is much larger than -|p at r=r b , and 
equation (A. 6) can be taken as 




Extending this solution to a general r 




or 




(A. 7a) 



and 



<j> = -n— 1 — S(t-r/a ) 
* 47Tr 00 



(A. 7b) 
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from equation (A. 2). Equations (A. 7) represent the pressure 
and velocity potential distributions of a stationary oscil- 
latory source of strength S. 



2. Moving Singularity 

Consider a source moving with uniform velocity U , with 
respect to the surrounding fluid, in the direction of the 
negative x axis. If Q(x,y,z,t) is the source distribution 
density, the continuity equation (2.1) becomes 

+ V-pU = Q 



and the linearized wave equation (A. 3) becomes 



V 2 p 



a^3t 2 



9Q 

3t 



(A. 8) 



This source distribution can be represented by 



Q(x,y,z,t) = p S(t) 6(x+U t) 6(y) 6(z) 



where 6 represents the familiar Dirac delta function. Making 
use of the linear relationship between p and $ in equation 
(A. 2), equation (A. 8) may then be expressed in terms of the 
velocity potential as 

2 - 

V 2 <f) - -kj — | = S(t) 6(x+U TO t) 6 (y ) 6 ( z ) (A. 9) 

at 9t 

A Lorentz transformation, scaled with respect to the compres- 
sibility factor 6, can now be used to reduce the right hand 
side of equation (A. 9) to that of a stationary source. The 
transformed coordinate system (prirr.ed) is defined as follows: 
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X' 



= ( x+U t ) 

b 



y' = y/B 



z’ = z/B 



(A. 10) 



u 



t' 8 i (t+-- X) 

e 4 

Performing this transformation, equation (A. 9) takes the 
form 



2-r 



V' 2 <t> ~4 ^4r = 4 S(t*)6(x')6(yM6(z*) 

B 
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(A. 11) 



Equation (A. 9) has therefore been reduced to the form of the 
wave equation for a stationary source in the primed coordi- 
nate system, which from equation (A. 7) has the solution 



<t>(r',t') = 44 ' S(f - r'/a w ) 



(A. 12) 



In terms of the space fixed coordinate system, the variables 
of equation (A. 12) are 






x+U t) 2 + e 2 (y 2 +z 2 ) 



u 



t' = 4 + -? x ) 

Transferring the solution (A. 12) to the wing fixed coordinate 
system now requires the further Galilean transformation 
x = x+U t, which produces the final result 



«><r,t) = ^ s[t +-±j (M^x-r)] 



(A. 13) 
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where 



r = ~\/x 2 + f^(y 2 + z 2 ) 

All variables are referenced to the wing fixed coordinate 
system. 

3. Wing Fixed Periodic Singularities 

If an harmonically oscillating source of unit strength 
(S = le^ wt ) is located at the origin of the wing fixed 
coordinate system, the velocity potential distribution given 
by equation (A. 13) is 

® s (x»y»z»t) = - exp jiu[t + — (M^x-r)]} (A.l4) 

a oo^ 

Considering now only the spatial variation of the velocity 
potential of a unit source located at a general point 
(x o ,y oJ z o ) (time dependence e lu ^ having been factored out) 

$ s (x,y ,z ) = ^ exp {i — ^ [m oo (x-x q ) - r] } (A. 15) 

a oo^ 

where 

R = ^(x-x Q ) 2 + 6 2 [(y-y 0 ) 2 + (z-z Q ) 2 ] 

The potential distribution of a dipole or doublet is 
related to that of a source by 

94 > 

*d = TrT 

where n is the direction of the axis of the dipole. There- 
fore from equation (A. 15), the velocity potential of a unit 
dipole oriented in the z direction is 



135 



3+, 

*d (x >5" 2) = W 



o 



-e 3 

^Jnk 



3 (z 'V 



1 + i- 



■7T2 R ) ex P( i 3[ M . (x - x o ) - R ]) 

a ~ e ' ae r& _u 



(A. 16) 



The downwash at a point on the wing (x,y,o) from such a 
doublet located in the plane of the wing ( x 0 >y o J°) is given 

by 



w(x,y,o) 
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o 



= 0 



or 
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1+ i — ^-pR expji — 
a co6 1 1 M 




(A. 17) 



where 



R = -\/(x-x o ) 2 + 6 2 (y-y Q ) 2 

This is the standard form of the kernel function in cartesian 
coordinates normally used in the velocity potential formula- 
tion of lifting surface theory. 



Kernel Function for Wing/Body Doublets 

To analyze the wing/body problem, a transformation from 
cartesian to cylindrical coordinates was made as discussed 
in Section IV. B and shown in Figure 8. Thus 

x = x 

y = r cos 0 
z = r sin 0 

relate the equations of the previous section to the (r,0,x) 
system. A unit doublet with axis in the z direction located 
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at ( r 0 >® 0 5 x 0 ) has a potential distribution, from equation 
(A.16), in the new coordinate system of 



4> d (r,e,x) = 



-8‘ 



^TT R 



^ (r sin 0-r Q sin 0 Q ) 1 1 + i 



to 
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exp j i — [m co (x-x q ) - r] } (A. 18) 



where 



a oo^ 



R = -\/(x-x ) 2 + 0 2 (r 2 +r 2 - 2rr cos <0-0 >) 
v o o o o 



For a doublet located on the wing at (r , 0 = 0 or Tr,x Q ), the 
downwash on the wing at (r,o,x) is given by 



u„ = — 



1 9<t> d 



0=0, 0 =0 or it 
’ o 



0 r 3 0 

which applied to equation (A.l8) results in 



U 0 = -^-3 

4ttR 3 



1 R 



exp |i — ^ [m (o (x-x q )- p]) (A. 19) 
& B 

OO 



where 



R = - n /(x-x q ) 2 + 0 2 (r^+r 2 + 2rr Q ) 



- for 0=0 
+ for 0 = TT 

This is of course the same as equation (A. 17) in cartesian 
coordinates except for R which indicates the coordinate 
trans formation . 



137 



The normalwash on a cylindrical body of radius r D , 

D 

coaxial with the x axis, at a control point, (r D ,0,x), due 

D 

to the wing doublet is given by 



u 



r 
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9r 
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0 =0 or it 
o 



or 




where 



(A. 20) 



R 



= "\/u- 
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(r D + r cos 0) 
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- for 0=0 
+ for 0 = it 

To determine the body singularity kernel function forms, 
the velocity potential of a source is expressed in cylindri- 
cal coordinates. From equation (A. 15) 

4> s (r , 0 , x ) = ^ exp (i [m co (x-x o ) - r] } (A. 21) 

1 a 8 
00 M 

where 

R = *\/(x-x o ) 2 + 8 2 (r' 1 +r 2 -2rr o cos <0-0 Q >) 
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The potential distribution of a doublet oriented in the 



radial direction and located at (r ,9 , x ) is obtained from 

o o o 

3 4 > 

v d 3r 

o 



which applied to equation (A. 21) gives 



<f>.(r,6,x) = - - fr-r cos (6-6 )"] fl + i - ? r] 
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(A. 22) 



with R as above. From this potential formulation, the 
following kernel functions can be obtained for a unit doublet 
with axis in the radial direction, located on the body sur- 
face at (r B , 9 q , x q ). 

The downwash on the wing at (r,o,x) is given by 



u = ilid 

u 6 r 36 



6=0, r =r 
o 
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or 




where 



R = -\/(x-x 0 ) 2 + 0 2 (r 2 +r B -2rr B cos 9 q ) 
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The normalwash on the body at (r fi ,0,x) is given by 
9<fc ' 
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(A. 2*1) 



where 



R 



= -^/(x-x ) 2 + 23 2 r 2 (1-cos <0-0 >) 

v O U O 



M = — r R ( 1-cos <0-0 >) 
o r p b o 

Equations (A.19)> (A. 20), (A. 23), (A. 2*1) represent the 
kernel function formulations used in the discrete potential 
element analysis of the wing and wing/body programs. 
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APPENDIX B 



INPUT INSTRUCTIONS FOR COMPUTER PROGRAMS 



Instructions for preparing input data for both the wing 
and wing/body programs are presented in this appendix. The 
field location and format for each input quantity is speci- 
fied. Any set of units may be used for geometric dimensions, 
displacements, and acoustic velocity as long as they are 
consistent. Any number of problems may be run in sequence, 
with the programs terminating when a new data set is not 
available to be read. 



1. Wing Program 

a. First card: FORMAT <3F10 . -4 ,2110 ,F10\ ) 

Column 1-10 11-20 21-30 4 0 50 51-60 

Name XM F A ICH JCH ALPH 

XM Mach number. 

F Circular frequency (rad/sec). 

A Speed of Sound. 

ICH Indicator for second card input data. 

= 0 for new data; program reads second card. 

= 1 for same data as previous problem; no 
second card required. 

JCH Indicator for wing displacement input data. 

= 0 for new data; program reads third and 

subsequent wing displacement cards. 

= 1 for same data as previous problem or 
for pitching motion; no wing displacement 
cards required. 

ALPH Amplitude of pure pitching motion. 

> 0 for pitching motion; program computes 
down wash for pitching mode (JCH must equal 1). 
= 0.0 for general wing motion; program com- 
putes downwash from wing displacement input. 



b. Second card: FORMAT( 415 ,5F10 . 4 ) 



Column 1-5 6-10 11-15 16-20 21-30 31-40 41-50 51-60 61-70 
Name N M MB NC XC XB BR Gil G22 

N Number of chordwise control points (maximum ten). 

M Number of spanwise control points (maximum ten). 

MB =0 for wing program. 

NC Significant wake length downstream from wing 

trailing edge (root chord lengths). 

XC Wing root chord. 

BR = 0 for wing program 

XB Wing semi-span. 

Gil Leading edge sweep angle (deg.). 

G22 Trailing edge sweep angle (deg.). 

Gil and G22 positive for downstream sweep. 



c. Third and subsequent wing displacement cards: 

FORMAT (7F10. 4) 

Column 1-10 11-20 21-30 31-40 41-50 51-60 61-70 

Name Z(l) Z(2) Z(3) Z(4) Z(5) Z (6) Z(7) 

Z ( 8) Z (9) Z(10) fourth card 

Z(ll) Z ( 12 ) Z (13) etc. fifth card 

Displacements (wing motion amplitudes) are read span- 
wise starting from the wing root at the leading edge. 
A new card must be started for each spanline of data 
(example above is for ten spanwise points). 



Wing/Body Program 

a. First card: same as wing program. 

b. Second card: FORMAT ( 415 » 3F10 ,15 ) 

Column 1-5 6-10 11-15 16-20 21-30 31-40 41-50 51-55 
Name N M MB NC XC XB BR MS 

N Number of wing chordwise control points 

(maximum ten). 

M Number of wing spanwise control points 

(maximum ten ) . 



MB 


Total number of control points 
quarter-circumference surface 


on body 
(maximum 


130) 


NC 


Significant wake length downstream from 
trailing edge (chord lengths). 


wing 


XC 


Wing chord. 






XB 


Wing semi-span. 






BR 


Body radius. 






MS 


Number of control points along 
circumference . 


body quarter 



Third and subsequent wing displacement cards: same 
as wing program. 



APPENDIX C 



COMPUTER PROGRAM FLOW DIAGRAMS 



START 




FIGURE 34 

WING PROGRAM - tAA\N FLOW DIAGRAM 
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FlGiURE 3S 

SUBROUTINE PRU3 FLOW OV/NoFAFA 
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FIGURE 3£> 

W\UG PROGRAVA -UlUCO FLOW DIGRAM 
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START 




FIGURE 37 



SUBROUTINE CO MAT" ELOW D\ AGRAM 




FIGURE 38 



SUBROUTINES SECUVA AND W\NGLM FLOW DIAGRAM 



Yes 




F\6URE 3S 

\a/)NG/SODY PROGRAM - YAA\N FLOW OlAGKAt^ 
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WlUCa/BODY PROGRAM - 'D\WCO ROW 0 \AGRAM 




F\GURE VvTAe 



SUbFOLTHNE UTHE.TA. FLOW 0\Av6RMA 
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START 




F\GURE 47_ 

SUBROUTIWE URAO.L FLOW Q\AGRAF^ 
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FIGURE 43 

SUBROUTINE UFAOZ. FLOW 0\AGRAF\ 
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APPENDIX E 



LISTING OF THE FORTRAN CODE FOR THE WING/BODY PROGRAM 
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